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ABSTRACT 

We investigate the dynamics of bow shock nebulae created by pulsars moving super¬ 
sonically through a partially ionized interstellar medium. A fraction of interstellar 
neutral hydrogen atoms penetrating into the tail region of a pulsar wind will undergo 
photo-ionization due to the UV light emitted by the nebula, with the resulting mass 
loading dramatically changing the flow dynamics of the light leptonic pulsar wind. Us¬ 
ing a quasi 1-D hydrodynamic model of both non-relativistic and relativistic flow, and 
focusing on scales much larger than the stand-off distance, we find that if a relatively 
small density of neutral hydrogen, as low as 10“^cm“^, penetrate inside the pulsar 
wind, this is sufficient to strongly affect the tail flow. Mass loading leads to the fast 
expansion of the pulsar wind tail, making the tail flow intrinsically non-stationary. 
The shapes predicted for the bow shock nebulae compare well with observations, both 
in Ha and X-rays. 
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1 INTRODUCTION 


It has been estimated that between 10% and 50% of 
pulsars are born with k i ck velocities Ujs 500 km s“ ^ 

llCordes fc ChernofilllQQ^ : lArzoumanian et al.l 200^ . These 
pulsars will escape from their associated supernova remnants 
into the coo ler, external interstellar m edium (ISM) in less 
than 20kyr (I Arzoumanian et al.ll2003) . As this time scale is 
sufficiently short, the pulsars are still capable of producing 
powerful relativistic winds. Furthermore, comparison with 
typical sound speeds in the ISM, Cs ,ism = 10 — 100kms”*^, 
shows that the pulsars are moving with highly supersonic 
velocities. The interaction of the pulsar’s wind with the 
ISM produces a bow shock nebula with an extended tail. 
If a pulsar is moving through a partially ionized medium, 
the bow shock nebula can be detected by the characteristic 
Ha emission resulting from the collisional and/or charge- 
exchange excitation of neutral hydrogen atoms in the post¬ 
shock flows and the subsequent e mission via bound-bound 
transitions dChevalier et alJIlQSd) . To date, nine such bow 
shock nebulae have been discovered, including three around 
7 -ray pulsars (iBrownsberger fc Romanill2014 ). 

Hydrodynami c (and hydromagnetic) models {e.g. 
lBucciantinill2002a ) of bow shock nebulae predict the forma¬ 
tion of a smooth two shock structure schematically shown 
in Fig. [U a forward shock in the ISM separated by a con¬ 


tact discontinuity from a termination shock in the pulsar 
wind. In the head of the nebula the shock and the contact 
discontinuity are situated at a distance given by Eq. ©, 
corresponding to the position where the ram pressure of the 
ISM balances the pulsar wind pressure. The flow structure 
in the head of the nebula is reasonably well understood, es¬ 
pecially in the limit of strong shocks {i.e., when the pulsar 
velocity is much larger than the ISM sound speed) and ne¬ 
glect ing the internal structure of the shocked layer llWilkinI 

Sfl 

These models further predict that the pulsar wind ter¬ 
minates at a Mach disk located approximately at a distance 
dback = Mns^o behind the pulsar (where Mns = Uns/cs.ism 
is the Mach number of the pulsar moving through the 
ISM). At approximately the same distance the oblique for¬ 
ward shock turns into a Mach cone with an opening angle 
~ 1/Mns- For distances larger than dback the flow in the 
tail of the nebula is smooth and nearly cylindrical, although 
some models predict the development of shear flow instabil¬ 
ities (e.g., Kelvin-Helmholtz instabilities). 

In contrast to these numerical models, Ha, radio and X- 
rays observations show that the morphologies of bow shock 
nebulae are significantly more complicated. More specifi¬ 
cally, observations reveal that the tails of bow shock nebulae 
have a highly irregular morphology, with Fig[2]showing four 
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Figure 1. Schematic illustration of a bow shock nebula prop¬ 
agating through fully ionized ISM, as seen in the rest frame of 
the pulsar. The dot-dashed rectangle shows the region zoomed in 

FigE] 


such examples. All these nebulae have a characteristic “head- 
and-shoulder” structure, with the smooth bow shock in the 
head not evolving into a quasi-conical or quasi-cylindrical 
shape, but instead showing a sudden sideways expansion(s). 
Arguably the most famous example is the Guitar nebula 
powered by the pulsar PSR B2224+65 (top left panel in 
Figi}. As the name suggests, this nebula has a guitar-like 
shape with a bright head, a faint neck, and a body consisting 
of several larger bubbles. 

Although the morphologies of these nebulae vary from 
source to source, there are a number of common features 
which, in our view, not only reflect the intrinsic dynamical 
properties of the flows, but which are also independent of 
the subtle details of both the pulsar winds {e.g., the relative 
orientation of the velocity and spin axis) and of the ISM. 
We stress the fact that all bow shock nebulae show qualita¬ 
tively similar morphological features not expected from sim¬ 
ple fluid models. In the X-ray and radio bands the tails show 
highly non-trivial morpho logies with quasi-periodic varia- 
tions in the intensity {e.g. iKargaltsev fc Pavlov!|200^ . For 
example, in the case of the Guitar nebula the tail shows 
quasi -periodic bubble-like structures (Ivan Kerkwiik fc Ingld 

l2008h . 

These peculiar tail shapes have be en interpreted as the 
result of density var iations in the ISM (iRomani et al.l[l997l : 
IVigelius et al.ll2007^ . However, based on the following con¬ 
siderations, we find this explanation unsatisfactory: (i) all 
tails show similar morphological variations (see Fig[2}; {ii) 
a common characteristic of these bow shock nebulae is that 
they are all highly symmetric with respect to the direction 
of motion of the pulsar - this is not expected if variations are 
due to the external medium; (Hi) morphological features in 
Ha, radio and X-rays are quasi-periodic - this is also not ex¬ 
pected from random ISM density variations. From these ob¬ 
servations we conclude that the peculiar morphological fea¬ 
tures result from the internal dynamics of the pulsar wind, 
rather than through inhomogeneities in the ISM. 


Ivan Kerkwiik fc Ingld ll2008h have also previously pro¬ 
posed that the morphology of the Guitar nebula could 
be explained by (unidentified) instabilities in the jet-like 
flow of pu lsar material away from the bow shock. Al¬ 
ternativ ely. iBucciantini fc Bandieral 1I2OOIII and iBucciantinil 
ll2002bll have suggested that the mass loading of pulsar wind 
nebulae may strongly affect their dynamics. These authors 
have shown that a non-negligible fraction of neutral atoms 
can cross the shocked ISM behind the bow shock without 
undergoing any interaction, thereby enabling these atoms 
to propagate into the pulsar wind region. Once inside the 
wind, neutral hydrogen can be ionized by UV or X pho¬ 
tons emitted by the nebula, and possibly by collisions with 
relativistic electrons and positrons, resulting in a net mass 
loading of the wind. 

_ In ord er to study this scen ario, IBucciantini fc Bandieral 

1I2OOIII and IBucciantinil ll2002bfl extended the thin -layer ap¬ 
proximation used to mo del cometary nebulae llBandieral 
I1993I : IWilkinIIl996l . I2OOOII . The thin-layer approximation is 
conceptually analogous to a 1-D model as it neglects the 
thickness of the nebula, while all quantities depend only on 
the distance from the apex. Despite the above-mentioned 
simplifications, these models provide a good description of 
the head region of the nebulae in terms of shape, hydrogen 
penetration length scale, and Ha luminosity, as was later 
confirmed by more accurate 2-D axisym metric simulations , 
both in the hydrody namic (HD) regime (lBucciantinill200^ : 
iGaensler et al.ll2004h an d in the relativistic ma gnetohydro- 
dynamic (MHD) regime I Bucciantini et ed]|2005ll . Using a 3- 
D model, IVigelius et al.l ( 2007ll was able to extend the study 
of these systems by also taking into account either a non- 
uniform ambient medium, or the anisotropy of the pulsar 
wind energy flux. However, none of these models are able to 
explain the peculiar morphology of the Ha emission often 
observed in the tail regions of bow shock nebulae. 

While the above-mentioned studies focused primarily 
on the head of bow shock nebulae, the aim of the present 
paper is to investigate the effect of neutral hydrogen on the 
tail region of these nebulae. The question we would like to 
investigate is whether the mass loading of neutral hydro¬ 
gen in the pulsar wind can explain the peculiar morphology 
observed at Ha, radio and X-ray energies. In order to fo¬ 
cus on the effect of mass loading on the evolution of bow 
shock nebulae, complications introduced by magnetic field 
pressure (and topology) are neglected in the present paper. 
These aspects are indeed necessary for a comprehensive and 
realistic treatment of the problem, and will be the subject 
of a future study. 

At this point the question arises as to whether one can 
use observations of the heliopshere to understand the prob¬ 
lem formulated in the previous paragraph. Although mass 
loading pl ays an important ro l e in the dynamics of the so - 
lar wind (iBaranov et al.l Il97ll : iBaranovI Il99(]l : IZankI Il999l l , 
there are a number of key differences between this scenario 
and the pulsar wind scenario. Firstly, the velocity of the 
Sun thr ough the ISM is, most likely, weakly sub-fast magne- 
tosonic (iMcComas et al.ll2012l l. whereas the pulsar’s motion 
is highly supersonic; secondly, the pulsar wind is very light 
- composed of lepton pairs - and one would therefore expect 
mass loading to have a greater effect; thirdly, most of the 
research related to mass loading in the solar wind concen- 
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Figure 2. Montage of Ha images of optical bow shocks associated 
with pulsar wind nebulae. Shown are J 2224+65, the so-called Gui¬ 
tar nebula dChatteriee &: Cordesll2002h. J07 42-2822. J2030+4415, 
and ,72124-3348 llBrownsbergerfc^omanill2014h . 


trated on flow in the head region, while we are interested in 
the large scale dynamics of the tail flows. 

The outline of the rest of the paper is as follows: in !j2] 
we discuss how neutral hydrogen penetrates into the pul¬ 
sar wind, summarising the possible interaction processes 
(charge-exchange, photo-ionization and collisions). In par¬ 
ticular, we show that the ionization of neutrals inside the 
wind is mainly due to UV photons emitted by the nebula. 
We stress here that in the whole paper with “neutrals” we 
will always refer to neutral hydrogen, even if, in principle, 
neutral Helium can play a similar role. We will also briefly 
argue at the end of § [5] why ions from the ISM should play 
no role in mass loading. The analytical, non-relativistic HD 
model used for the study is presented in ^ and Jjd] while 
in 3S]we develop a similar model but for a relativistic wind. 
In !jB]we provide a simple visual model fit to the bow shock 
nebula PSR J0742-2822, while we qualitatively discuss in m 
how the presence of a magnetic field inside the wind can 
modify the wind dynamics. Lastly, a summary of the main 
results can be found in SJH] 


2 STRUCTURE AND RELEVANT SCALES OF 
A PWNe CONFINED BY A PARTIALLY 
IONIZED MEDIUM 

Fig.fT]) shows the typical structure produced by a pulsar 
(or a normal star) moving supersonically through the ISM, 
as seen in the rest frame of the pulsar, when the effect of 
neutrals is not taken into account. Such bow shock struc¬ 
tures are preferentially produced when the pulsar propa¬ 
gates through the warm phase of the ISM, whose typical 
temperature is T ~ 6 • 10^-10^ K. In this case the sound 
speed is Cs ~ 10 km/s, hence the pulsar’s Mach number is 
';$> 1. Conversely, the hot phase of the ISM has T ~ 10® 
K and Cs ~ 100 km/s, implying a Mach number close to 
one. In the hot phase the ISM is totally ionised, while in 
the worm one the presence of neutral Hydrogen cannot be 


neglected. In fact the typical density of the warm ISM is 
0.2-0.5 cm“® and the ionized fraction is estimated to range 
between 0.007-0.05 for the warm neutral medium (WNM), 
and 0.6-0.9 for th e ionised neutral medium (WIM) (see, e.g., 
Ijean et~^ I I2OO9I. Table 1). 

The distance, do, between the pulsar and the contact 
discontinuity (CD) (formed between the shocked ISM and 
the shocked wind) is obtained by equating the wind pressure 
with the bulk pressure of the ISM: 


do = 




(dvrUjIgPisMc/ 




o ,, inl6 /'1/2 T7-1 -1/2 

3 X 10 34 V 300 lliSM,-! '"'21 , 


( 1 ) 

where T™ = erg s ^ is the pulsar luminosity, 

Uns = SOOVsookms”^ is its peculiar velocity, and pisM = 
mpTiisM is the density of the dragged component of the ambi¬ 
ent medium, expressed in units of nisM = 0.1 nisM.-i cm“®. 

In order for mass loading to play a role in the dynamic 
evolution of bow shock nebulae, neutrals are required to 
cross A (the distance between the bow shock and the contact 
discontinuity) and penetrate into the wind region. Hence the 
interacti on length in s ide th e shocked ISM, A, must be larger 
than A. IChen et al.l (Il996l l estim ated that A = 5/1 6 dp, a 
value that was later confirmed bv iBucciantini et al.l ll2005l i 
using numerical simulations (note that A is smaller than 
do). 


At this stage it is important to note that different phys¬ 
ical processes will lead to the interaction of neutrals in the 
shocked ISM than in the pulsar wind. In the next section 
it will be shown that a significant amount of neutrals can 
undergo charge exchange (CE) with protons in the shocked 
ISM and collisional ionization with electrons (specifically in 
the head of the bow shock system), while the collisional 
ionization of neutrals by protons can be neglected when 
Vns JS 300kms“^. Despite the neutrals undergoing CE and 
ionization, it will further be shown that a dynamically im¬ 
portant fraction of neutrals can still penetrate into the pul¬ 
sar wind in their original state. 

The pulsar wind most likely consists of o nly electrons 
and positrons lIRuderman fc Sutherland! 1 1975l l. and CE is 
therefore not possible. Rather, ionization of neutrals can 
only occur through photo-ionization or through the colli¬ 
sion with relativistic electrons and positrons. The former 
process is discussed in 112.31 where it is shown that a signif¬ 
icant amount of neutrals can be photo-ionized through the 
non-thermal emission emitted by the pulsar wind, while in 
112.4l we show that collisional ionization can be neglected. As 
photo-ionization is important for the pulsar wind, one may 
ask whether this process is also important for the shocked 
ISM in the head of the nebula. Section l)2.2l will therefore be 
used to discuss when photo-ionization in the shocked ISM 
can be neglected. Lastly, for the sake of completeness, we 
will discuss in 112.5l whv the ionization of neutrals inside the 
unshocked pulsar wind can be neglected. 

One may wonder whether ions from the ISM, rather 
than neutrals, can penetrate inside the wind directly, re¬ 
sulting in a net mass loading of the wind. The main is¬ 
sue of this hypothesis is that ions are attached to magnetic 
field lines which does not cross the contact discontinuity. In 
principle ions could diffuse perpendicular to the held lines 
and enter inside the wind. Nevertheless a simple estimate 
excludes this possibility: using the Bohm diffusion coeffl- 
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cient, Db = wl/S, as an upper limit for the perpendicular 
diffusion, we can estimate the time needed for a thermal 
proton to cross the typical distance do, which is given by 
tcross = do/Ds- Using B = 1/rG, do = 10^® cm and v = 100 
km s~^, we get tcross ~ 10® yr. Hence the diffusion of ions is 
too slow and can be neglected. An alternative way for ions to 
penetrate inside the wind is through shear flow instabilities 
that can develop at the contact discontinuity. Nevertheless, 
the effects of instabilities will not be limited to the simple 
injection of ions, but will mix the ISM materials with the 
pulsar wind resulting in a complex tale structure and prob¬ 
ably to its disruption. In our knowledge this scenario has 
never been studied in details, and it is far beyond the aim 
of the present paper. 


2.1 Interaction of neutrals in the shocked ISM 


After crossing the bow shock, neutral atoms can interact 
with the shocked protons. The interaction length is given by 


_Uns_ 

A^ionUisMr c (rr (nrel)ni.el) 


( 2 ) 


where Aion is the ionization fraction of the ISM, is the 
compression ratio of the bow shock, cr(iirei) is the rele¬ 
vant cross section of the process under consideration, and 
(cVrei) is the collision rate averaged over the ion distribu¬ 
tion function. When the ion distribution is a Maxwellian, 
(aVrei) is well approximated (within 20%) by the expression 
( Zank et al.lll99a : iBlasi et al.ll2012h 

(aVrel) ~ a{U*)U* , (3) 


where 


U* 


8 2kBT 


(4) 


is the average, relative speed between the incoming hydro¬ 
gen atom and ions (T is the temperature of the shocked 
ISM determined by the Rankine-Hugoniot jump conditions, 
assumed to be ':$> Tism). Using the fiducial values nisM = 
0.1cm“® and Uns = 300kms“^, together with an ioniza¬ 
tion fraction of 90% and rc = 4 (the typical value for strong 
shocks), leads to the following estimates for the mean free 
paths 


‘ion,p ^ 

s 3.0 X 10®° cm 

(5) 

^ion,e ^ 

s 2.2 X 10^° cm 

(6) 

Ace " 

s 1.5 X 10^® cm 

(7) 


for the ionization due to collisions with protons, electrons 
and CE, respectively. Note that Aion.e has been calculated 
under the assumption that the electrons downstream of the 
shock equilibrate rapidly with protons, thereby acquiring 
the same temperature. If this assumption does not hold, the 
collisional length scale for ionization due to electrons can 
become much larger than the value reported in Eq. ((6|. In 
addition, the values are to be taken as lower limits 

as they are valid just ahead of the nebula, where the com¬ 
pression ratio and the temperature obtain their maximum 
values. 

From these estimates it follows that only a negligible 
fraction of the neutral hydrogen will be collisionally ionized 
by electrons, whereas a significant fraction of neutrals will 


undergo CE. The neutrals resulting from a CE event will 
have a bulk speed and a temperature that are close to that 
of the protons in the shocked ISM. This implies that the 
newly formed neutrals tend to be dragged with the shocked 
protons along a direction parallel to the contact disconti¬ 
nuity. Nevertheless, the CE process produces a diffusion of 
neutrals in the nose of the nebula and it may still be possi¬ 
ble for the newly formed neutrals to enter the wind region, 
provided that their diffusion velocity perpendicular to the 
contact discontinuity is of the same order or larger than 
their velocity parallel to the contact discontinuity. This is a 
complication that will not be addressed in the present paper 
but is essential to estimate the correct amount of neutrals 
that can penetrate into the wind. 

Although a large number of neutrals will be lost due 
to CE, a minimum fraction of neutrals proportional to 
exp[— A/(Ace -|- Aion)] will cross the shocked ISM region 
without suffering any interaction and will enter the pulsar 
wind in their original state. These neutrals will not influence 
the wind structure until they are ionized, either through col¬ 
lisions with relativistic electrons and positrons, or through 
photo-ionization with photons emitted by the nebula or by 
the pulsar. In l)2.3l we show that photo-ionization is the dom¬ 
inant process in pulsar wind nebulae, while we demonstrate 
in 112.4l that collisional ionization can be neglected. However, 
we first discuss the photo-ionization ahead of the nebula and 
in the shocked ISM. 


2.2 Photo-ionization outside the pulsar wind 

There are three different sources of photons to account for: 
thermal and non-thermal radiation from the pulsar, and 
non-thermal radiation from the nebula. The thermal ra- 
diati on from the pulsar has been shown to be negligible 
Isee iBucciantini fc Bandieral 120011 . Eq.(25) and discussion 
below). On the other hand, non-thermal emission from both 
the pulsar and the nebula can play a role as the non-thermal 
pulsar luminosity is generally comparable to the luminosity 
of the nebula. 

When the radiation emitted by the nebula has a suf¬ 
ficient luminosity, incoming neutrals can be ionized before 
crossing the bow shock, and one would consequently not ex¬ 
pect any effect from mass loading. Conversely, the require¬ 
ment that hydrogen atoms are not fully ionized at the bow 
shock imposes an upper limit on the total ionizing luminos¬ 
ity. We d erive this limit closely following a similar derivation 
given by Ivan Kerkwiik fc Kulkarnil ll200llf . 

Assuming a spherical symmetry for the emission emit¬ 
ted by the nebula, the ionization fraction at a distance r 
from the center of the nebula is 

di+{r) iVphe-"W , 

— (1 €+)<rph 5+nearec , (8) 

where Aph is the number of ionizing photons emitted by the 
nebula per unit time, and (tph = N~^ ap\^{v)nph{v)dv is 
the photo-ionization cross section averaged over the photon 
distribution. The photo-ionization cross section is 

cTphG) = 64a'^^aT {I/huY^'^ = 10“^®(/izz/Ryd)“^^® cm®, 

(9) 

with (Tt the Thompson cross section and I — 1 Ryd the 
ionization potential. As Cph decreases rapidly with photon 
energy, the only relevant photons are those with energy > I. 
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J ulsar, non- 

1 I 2 OI 4 I) . ( 2 ) 

lAbdo et alj ll2013l') . (3) iHui &: Becked ll2007l) , (4) I Romani et alj ll2010h . (5) IStappers et alj ll2003h . (6) iHui fc Becker! ll2006h . The values 
of Lx marked with an * are calculated using the relation log^QLx.pwn = l-^Hog^^QL — 21.4 derived by iKargaltsev e^^ I ll2012h from 
Chandra observations. 


Pulsar 

E 34 

-^X,pul 

Ax,pwn 

rX,pwn 

msM 

Reference 


[ergs-l] 

[ergs-l] 

[ergs-l] 


[cm~®] 


J0437-4715 

0.55 

2.4 X 10®° 

3.5 X 10®° * 

- 

0.21 

1,2 

J0742-2822 

19.0 

< 9.6 X 10®° 

7.4 X 10®i * 

- 

0.28 

1,2 

J1509-5850 

68.2 

2.4 X 10®® 

(1.3 - 2.7) X 10®® 

-1 o-)-U.» 
-‘-•'^-0.4 

6.14 

1 , 2,3 

J1741-2054 

12.6 

3.5 X 10®° 

4.0 X 10®i * 

1.6±0.2 

1.44 

1 , 2 , 4 

J1856-3754 

3 X 10"^ 

- 

4.2 X 10®"^ * 

- 

0.05 

1, 2 

J1959-I-2048 

21.9 

5.2 X 10®i 

1.9 X 10®i 

1.5 ±0.5 

0.02 

1, 2, 5 

J2030-I-4415 

2.90 

2.7 X 10®i 

4.3 X 10®° * 

- 

2.69 

1, 2 

J2124-3358 

0.68 

8.6 X 10®° 

10®° 

2.2 ±0.4 

0.47 

1, 2, 6 

J2225-I-6535 

0.16 

- 

5.5 X 10®® * 

- 

1.43 

1, 2 


Solving Eq.® in the rest frame of the pulsar where the 
plasma is moving in the positive direction along the ®-axis 
with velocity Vns, one may write dx = Vt^sdt. If the ex¬ 
tinction and recombination of protons are neglected, Eq.®, 
written in cylindrical coordinates {x,p), simplifies to 

d^-^- UphfVph dx finf 

1 — dvrCNS + x"^ 

The solution is straightforward and reads 


TT 

- arctan 

2 

( 11 ) 

where is the original ionization fraction of the ISM 
far from the nebula and ro = frphVph/(47rI/Ns) defines the 
typical distance where ionization is ef fective. As pointed 
out by Ivan Kerkwiik fc Kulkarnil (1200 ill , Eq. (fTTIl describes 
an ionization fraction which, for fixed impact parameter p, 
smoothly increases as one approaches the nebula from the 
left in Fig. 0 , strongly increases for p ~ ro, before approach¬ 
ing the asymptotic value 1 — (1 — as one moves 

away from the pulsar towards the right. 

In order to estimate the value of ro, it is useful to express 
the product (tph • Vph using the X-ray luminosity of the neb¬ 
ula, Cx, which we define as the luminosity in the Chandra 
energy range, *.e., between ei = 0.5 keV and £2 = 8 keV. 
Assuming that the non-thermal emission from the nebula 
ranging between UV and X-ray energies is a single power 
law with a photon index E, the value of ro can be written as 


^+{x) = 1 - (1 - C+,o)exp 



ro = = 1.33 • 10^®J-(r)£x.3oE3oo cm, (12) 

47r\/NS 

where Cx,30 is the X-ray luminosity in units of 10®® ergs~^ 
and V 300 is the pulsar speed in units of 300kms“^, while J- 
is a dimensionless function that depends only on the photon 
index 


j^(r) 


j9/2 


2T-4 [e 


- 5 / 2 - 




2r -b 5 [£2-r]' 


(13) 


The typical values of T inferre d for PWNe detected in X - 
rays are T ^ l.iQ (see also, e. 0 . IKargaltsev fc Pavlovll2008ll . 


^ It should be kept in mind that a value of F ~ 1.5 is expected if 


It follows that T ^ 2 y, 10“®, with T reaching unity when 
r = 2.55. Ea. (ll2ll shows that ro can easily be smaller than 
the stand-off distance, do, implying that the majority of neu¬ 
trals reach the bow shock without having been ionized. It 
is useful to express the condition ro < do in terms of an 
upper limit for the X-ray efficiency of the nebula, defined as 
px = CxICw'- 

ro<do ^ Vx< 10-" .F(r)-i Vsm'-i • (14) 

Typical measured values for PWNe are px ~ 10“®, although 
large deviations from this value are observed. For Ha emit¬ 
ting bow shock nebulae it is known that hydrogen atoms 
are present at the position of the bow shock, therefore the 
inequality (US has to be satisfied. In fact, for PWNe emit¬ 
ting both Ha bow shock and X-ray radiation, the estimated 
value of px is closer to 10“^ (see values reported in Table®, 
and the inequality m is thus satisfied. 


2.3 Photo-ionization inside the pulsar wind 


The photo-ionization of atoms inside the wind determines 
the rate of mass loading. In order to estimate this rate, the 
value of the photon density inside the nebula is required. 
For simplicity it is assumed that the wind is a cylinder with 
a cross section Tidp and a length R, emitting radiation uni¬ 
formly. The photon number density inside the wind is then 
given by 


^ph,w 


Vph(0 

rrd'^Rc ’ 


(15) 


The value of R can vary greatly from one bow shock nebula 
to another, but in general i? 3> do- The quantity (1) is the 
mean length traversed by a photon before escaping the neb¬ 
ula. When the presence of neutrals is neglected, PWNe are 
transparent to UV radiation because the le~ pair den¬ 
sity is very small and the Compton scattering mean free 
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the particles responsible for the non-thermal X-rays are produced 
through Fermi acceleration. 
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path is greater than the size of the nebul£0- Hence for geo¬ 
metrical reasons we can assnme (Z) ~ do. If we account for 
the presence of neutrals inside the wind, the mean free path 
of photons due to ionization is Zmfp = l/(crphnjv), where 
(jph is given by Eq.(l9]) and njv is the neutral density inside 
the wind. We do expect njv < njv.iSM ~ 0.1 cm“®, hence 
Zmfp > 10^^ cm. In conclusion, also when we account for 
the presence of neutrals, the assumption (Z) ~ do remains a 
good estimate and the ratio R/{1) is S> 1. Using Eq. (I15II . 
the ionization length of neutrals inside the wind is estimated 
to be 

, Uns 

^ph — _ 

^ph ,w <^phC 

= 3.2 • 10'® 14oo cm, (16) 

where the second equality has been obtained using Eqs. © 
and Ga. One important comment is in order: the ionization 
length scale, Aph, is ultimately not the length scale that 
determines whether neutrals will influence the dynamics of 
the wind. Aph is an estimate of the length scale required to 
ionize the largest fraction of neutrals. The quantity more 
pertinent to the dynamics of the wind is the relative change 
in the density of the wind induced by mass loading. Since 
the pulsar wind is light, consisting of electron-positron pairs, 
only a small fraction of neutrals are required to be ionized in 
order for these neutrals to become dynamically important. 
The more important parameter is therefore the length scale 
where the dynamics of the wind is dominated by the ionized 
protons rather than by the electron-positron pairs. Thus, we 
define AML.ph as the length scale where the loaded mass is 
equal to the initial mass of the wind, i.e., 


A 


ML,ph 


Pe Wvind 
PN ^ph,wf^ph C 


TheTTle V^ind 
UNTUp VnS 


Aph . 


(17) 


A numerical estimate of Aml.pH can be found using Eqs. m 
to evaluate Aph, and estimating the electron density of the 
wind. Tie, by equating the luminosity of the pulsar with the 
energy flux in the shocked wind, i.e., M = CwIGc?) = 
nemeirdoC. The result reads 


AhiL.ph = 1.32 • 10 


16 Uvind R / PX 

~lj) Vl(P 


)- 


n-^ 


(18) 


It is seen that AiviL.ph can be much smaller than do due to 
the fact that its value is inversely proportional to the elec¬ 
tron Lorentz factor, whose typical value is 7 « 10® — 10®. 
Even for Lorentz factor as low as 10®, a v alue estimated from 
the modelling of some pulsars (see, e.g.. lDubus et af]l2015ll . 
AML.ph remains smaller than cZq. It will be shown in 33] that 
Ea. (ll7|) is the correct definition of the length scale that de¬ 
termines whether mass loading will influence the dynamics 
of a non-relativistic wind. On the other hand, for a rela¬ 
tivistic wind the definition of Aml has to be modified as it 
becomes necessary to take into account the relativistic iner¬ 
tia of the wind. We will discuss this issue in 33]where we will 
derive the relativistic version of Eq. 113 using a relativistic 


® This conclusion could be different only for very young PWNe 
which have a larger density (see lAtovan fc AharoniarJ Il996l . 
for details). 


HD model. Moreover, we will show that there is no signifi¬ 
cant difference between the relativistic and non-relativistic 
results, apart from the different definition of Aml. 


2.4 Collisional ionization inside the pnlsar wind 

The collisiona l ionization lengt h scale depends on the Bethe 
cross section (iKim et al.ll2000l) 

aBethe(7) = [M^ (ln(7"/3") - + Cr] , (19) 

where 7 is the Lorentz factor of the electrons (and positrons) 
and P = {1 — 7 ~^)'^®. The two constants and Cr are 
related to the atomic form factors of the target, and are 
independent of the incoming particle ene rgy. For hydroge n 
atoms = 0.3013 and Cr = 4.3019 (iKim et alJl 20 ()ol ). 
Using the expression above, together with the typical value 
7 = 7510 ®, furnishes the collisional ionization mean free path 
in the pulsar wind: 

Acoi =- ' = 3.86X 10®® Vioo nisM,-i 75 cm . (20) 

As discussed in the previous section, the electron density 
of the wind is estimated by equating the luminosity of the 
pulsar with the energy flux in the shocked wind, i.e., M = 
CwIGc?) = UemeT^doC (thereby ensuring that the dynamics 
of the pulsar wind is not dominated by the magnetic field). 
Using the same values for the parameters as those used in 
Eq.GJ leads to the estimate Ue ~ 8 -10 ® 75 ' uism.-i cm ®. 
Note that in the relativistic regime (7 ^ 1) (TBethe increases 
only logarithmically with energy, and that increasing 7 from 
10® to 10® implies a decrease in Acoi of only ~ 40%. 

Similar to the case of photo-ionization, (see the discus¬ 
sion after Eg. (1161) 1 we define Aml, col as the typical mass¬ 
loading scale where the loaded mass density is equal to the 
initial mass density of the wind, i.e., 

. Pe Uwind f. — i ol9 Uvind I 

Aml.coI = --= 2.7 X lO-cm. 

PJV neO-Bethe(7) C C nN,-4 

( 21 ) 

It is again stressed that the above expression is for a non- 
relativistic wind and that the correct expression for the en¬ 
thalpy should be used for a relativistic wind. Comparison 
of Ea. (l2ll) and Ea. (ll 8 l) shows that photo-ionization is con¬ 
siderably more effective than collisional ionization, even for 
very low luminosity nebulae. Collisional ionization can thus 
be safely neglected in all realistic cases. 


2.5 Ionization inside the free wind region 

There is a finite probability that neutrals can be ionized 
inside the region of the unshocked wind (shaded region in 
Fig-©)- These ions could change the wind dynamics but, 
as we argue below, this process can be neglected. In fact, for 
parameters typical to bow shock nebulae, the Larmor radius 
of these ions, calculated in the rest frame of the wind, is 
comparable or larger than the size of the unshocked wind 
region. 

To show that this is, indeed, the case, an expression re¬ 
lating the termination shock radius and the wind luminosity 
is required. Assuming that the energy fluxes of the magnetic 
field and particles are similar, the wind luminosity can be 
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Figure 3. Sketch of the simplified model used to study the effect 
of mass loading in the tail of the bow shock nebula. This Figure 
represents a zoom in of the dot-dashed rectangle of Fig. Q. Note 
that the presence of the bow shock is neglected. The dash-dotted 
arrows show the regions where the neutrals in scenario B pen¬ 
etrate into the wind from the side of the contact discontinuity. 
This is in contrast to scenario A where neutrals only penetrate 
into the wind in a cylindrical region with cross section Ai. 

estimated as where Rt is the radius of the 

termination shock and B is the magnetic field strength at 
the termination shock. In the rest frame of the wind the 
newly formed ion have a bulk Lorentz factor equal to the 
wind Lorentz factor, 7 u,, hence the Larmor radius is: 

_ 'ym'iripP _ 7™mp 

eB - ^ 


For the sake of simplicity the presence of the bow shock is 
also neglected, but note that possible effects of mass load¬ 
ing on the shape of the bow shock will be discussed in 
It is further assumed that the external medium has a spa¬ 
tially independent velocity, Vo, pressure, Pq, and density, po- 
Lastly, it is assumed that both the ISM and the wind are 
non-relativistic with an adiabatic index 7 g = 5/3. In the 
next section we will develop the model for a wind with a 
relativistic temperature moving with a non relativistic bulk 
speed. The latter case is the more realistic scenario for a 
pulsar wind, whereas the former case applies more to stel¬ 
lar wind^. We anticipate that the main difference between 
these two cases lies solely in the different values used for the 
enthalpy. While this changes the dynamical length-scale sig¬ 
nificantly, it does not change the qualitative behavior of the 
solution. Finally, note that the steady-state nature of the 
system is not guaranteed. In fact, we will show that scenar¬ 
ios exist where the steady-state assumption is most likely 
violated. 

The conservation equations for mass, momentum and 
energy for a quasi 1-D system, written in the rest frame of 
the pulsar, are 



dj, [puA] = qA ', 

(24) 


9a: [pu^A\ + Ad^^t 

^ = qA'Vo , 

(25) 

dx 


= qA'vS/2. 

(26) 


Note that p = pe + Pp is the total density of the wind and 
that the quantities on the right-hand side of Eqs. (IMI-dSl 
represent the mass, momentum and energy flux due to the 
ionization of neutrals, respectively. The mass loading term 
is given by 


and the ratio between tl and the size of the free wind region 
can be expressed as 


Tl 

Rt 




5 76 -C 


- 1/2 
u;,34 • 


(23) 


The ionization of neutrals in the unshocked wind can pro¬ 
duce ions with tl > Rt that will escape the free wind re¬ 
gion, suffering only a small deflection. The requirement for 
the production of ions in the free wind regions is a high bulk 
Lorentz factor and a low spin-down power Cw Consequently 
we neglect the role of these ions and we will only account 
for neutrals ionized in the shocked wind. 


3 HYDRODYNAMIC MODEL FOR 
NON-RELATIVISTIC WIND 

To study the effect of mass loading in the tails of bow shock 
nebulae, the steady-state conservation equations for mass, 
momentum and energy are solved in a quasi 1-D approxi¬ 
mation. Fig. m shows a sketch of the model and represents 
an idealization of the region enclose in the dot-dashed rect¬ 
angle of Fig. ©■ The quasi 1-D approximation implies that 
the transverse cross section. A, of the flow can change, but 
that all the characteristic quantities of the wind, i.e., the 
velocity, u, density, p and pressure, P, are assumed to be a 
function of the position x only. This approach further implies 
that any internal structures are neglected, in particular the 
free wind region and the termination shock shown in Fig. 0. 


q = n{me + rup) . (27) 

As we showed in 1)2.31 n is determined predominantly by 
photo-ionization, allowing one to write 

h = UNUpuaptiC, (28) 

where njv is the local density of neutrals, Uph is the pho¬ 
ton density as seen in the rest frame of neutrals, and (fph 
is the photo-ionization cross section averaged over the pho¬ 
ton distribution. For the sake of simplicity we assume that 
both njv and Uph are constant along x. The spatial indepen¬ 
dence of riN is a good approximation when x Aph {i.e., 
the photo-ionization length scale given in Ea. (ll6ll '). while it 
becomes necessary to include the spatial dependence of njv 
when X > Xph. This can be done with a simple change of 
variables as we show in Appendix Conversely, a correct 
evaluation of the spatial dependence of the photon density is 
non-trivial and requires a detailed analysis of the specific ob¬ 
ject one wants to study. Here we avoid such a complication 
by assuming a constant value. 

It is also assumed that the incoming neutrals have the 
same temperature as the ISM (« 10"^ K), and they can thus 
be considered as being cold with respect to the wind in the 
nebula. As a result of this assumption, the momentum and 


For stellar winds it is necessary to also take into account the 
charge-exchange that occurs between neutral hydrogen coming 
from the ISM and protons inside the wind. 
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energy injection terms on the right-hand side of Eg. (1251) and 
(1261) only contain the contribution due to the bulk motion. 

The variable A' represents the effective area crossed by 
neutrals. In ^it was noted that neutrals can penetrate into 
the wind by crossing the head of the nebula. It is thus re¬ 
quired that the effective area should be A' ^ Ai = Txd^. 
However, in reality the situation is more complicated. We 
also noted that the CE process produces a diffusion of neu¬ 
trals in the shocked ISM and that this can lead to neutrals 
penetrating into the wind from the side of the contact dis¬ 
continuity. A rough estimate of this process, considering a 
constant diffusion coefficient, would result in a mass load¬ 
ing that is proportional to A{xY^'^. Moreover, it will subse¬ 
quently be shown that the effect of mass loading is to expand 
the cross section of the wind. When this happens the dis¬ 
tance between the bow shock and the contact discontinuity 
is reduced, thereby further increasing the probability of neu¬ 
trals penetrating from the side of the contact discontinuity, 
specifically in the tail region of the pulsar wind. This re¬ 
sults in a mass loading that is proportional to the total area 
A(x). A realistic solution therefore requires a description of 
the bow shock geometry and how the neutrals interact in 
the shocked ISM. Such a complication necessitates the use 
of a 2-D simulation, and is beyond the scope of the present 
work. In this study the two opposite situations, A' = A\ 
(scenario A) and A' oc A (scenario B), are investigated as 
we expect a realistic situation to be bracketed between these 
two scenarios. 

A consequence of the steady-state assumption is that 
the pressure is required to be constant everywhere, i.e., P = 
Pq and dxP = 0, and it is thus possible to simplify Eqs. (|25l) 
and (1261) . 

As a first step to finding the solution, Eq. (1241) can be 
substituted into the right-hand side of Eqs. (1251) and (1261) . 
leading to 


-puA [u^ -Vq) + 


dx [puA{u — Vo)] = 0 

^gPouA 


Ig - 1 


= 0 . 


(29) 

(30) 


These equations define two constants of the system which 
can be used to write down two expressions, one for p and the 
other for A, as functions of the velocity u only. Evaluating 
the constants of the system at the position x = xi = 0, 
where the boundary conditions are defined as w = Mi, 
A = Ai and p — pi, the following expression for the area is 
obtained: 

^ ^ , (m - u){ui - Vo)Gg - l)Mf 
Ai u 2uui ’ ^ ^ 

and the following expression for the density: 


u-Vo {ui-u){u-Vo){jg-l)MG \ " 

ui-Vo'^ 2ul \ ’ 

(32) 

where Mi = ui/csi is the initial Mach number of the wind 
and Cs ~ y/^gP!p is the sound speed. Note that all quan¬ 
tities with the subscript i refer to values at the boundary 


p(m) 

Pi 


Xl. 

The next step requires finding an expression for the ve¬ 
locity. Differentiating Eq. dSSl) by parts, and using Eq. (l24l) 
leads to 


dxU = —q 


A'{u-Vo) 

puA 


(33) 


This differential equation can easily be integrated by parts 
from a:i = 0 to x, after the quantity puA has been sub¬ 
stituted using Eg. (1291) . The solution for scenario A, where 
A' = Ai, is straightforward and reads: 


u{x) 


Xui + xVq 

X 


(34) 


where we have introduced the length scale A = uipi/q. If 
q is calculated using the photo-ionization cross section, A 
corresponds exactly to the definition provided in Eq. (O- 
An important and noteworthy result is that the solution 
(Oil) does not depend on the initial Mach number of the 
wind. Eq. (Oil) further shows that u(x) decreases monoton- 
ically from m(0) = ui to m(oo) = Vo, and that the typical 
length scale for this transition is Aui/Vb. The solution for 
scenario B {i.e., A' = A) is slightly more complicated but 
can be expressed in an implicit form as follows: 


X Ml — M 


-I- F 2 In 


Ml — Vo 2u\ + a{ui 
u — Vo 2uf 


m) 


(35) 


where 


El = 


2mi Vq 


F 2 = 


2u\ + a(Mi — Vo) 
2mi(mi — Vb)(2Mi -I- a) 
[2u'l -\- a{ui — Vo)f 


and 


(36) 

(37) 


a = Ml{ui-Vo){^g-l). (38) 

The solution (1351) also shows that m(oo) = Vq. This can 
be readily understood: when the amount of mass loaded is 
much larger then the initial mass of the wind, the final state 
of the wind will be the same as the initial state of the neu¬ 
trals. This is confirmed by the behaviour of Cs, which is also 
a monotonic decreasing function of 2 ;, as well as by the result 
u ^ Vo ^ Cs(m) —>■ 0 , which one expects from the assump¬ 
tion of cold neutrals. The continual process of mass loading 
therefore decelerates the wind from mi to Vq. The main dif¬ 
ference between solutions (1341) and (1351) is the typical length 
scale where this transition occurs: in the former case it is 
Xa ~ Ami/Vo, whereas in the latter case it is As ~ A. This 
result once again has a simple explanation: the amount of 
mass loaded in scenario A is a factor Ai/A smaller compared 
to the amount of mass loaded in scenario B. For x ^ 00 the 
solution (EH shows that the cross section of the flow in¬ 
creases asymptotically to the value 

lim A{u) = Ai — — Mi\ , (39) 

X—^oo VO \ 4 j 

where this asymptotic value is the same for both scenarios 
A and B. As we expect the transition length scale to be 
proportional to the loaded mass, we have A a — XbAo^/Ai ~ 
XbUi/Vo. 

The solutions for u{x), M{x), and A{x) corresponding 
to scenario A are shown in Fig. (jd)), and the solutions corre¬ 
sponding to scenario B in Fig. ©• For the latter case three 
panels are presented showing the results for three different 
values of the Mach number, while for scenario A only the 
case Ml = 1 is shown as the solution depends only weakly 
on the value of Mi. All these solutions are obtained by using 
the benchmark values summarised in Table [2] These values 
lead to an asymptotic expansion of A^o/Ai ~ 1000, which, 
in turn, translates into a radial expansion of a factor ~ 30. 
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Figure 4. Structure of the wind when mass loading occurs in 
scenario A, as calculated using the non-relativistic model. The 
plot shows, as a function of the position, the wind velocity divided 
by Vq (solid line), the Mach number, M (dashed thin line), and 
the expansion velocity Ur normalised to Uri, as given by Eq. lf44ll 
(dot-dashed line). The plot also shows the normalised area, A/A\ 
(dotted line) and the ISM sound speed normalized to Vq (thin 
dot-dashed line). The initial wind velocity is ui/Vq = 1000 and 
the initial Mach number of the wind is Mi = 1. Note that the 
spatial coordinate, x, is normalised to the mass loading length 
scale A. 


Comparison of Figs. and ©) confirms that the tran¬ 
sition in scenario B occurs much faster than in scenario A. 
However, apart from the different scale lengths of the tran¬ 
sition, the two scenarios are very similar. Fig. further 
shows that when the initial Mach number increases from 0.3 
to 3, the transition occurs faster by a factor ~ 3.5. 

A salient feature of Figs. and ® is that there is no 
qualitative difference between the solutions of a supersonic 
and a subsonic wind. This is a remarkable result given that 
the state of the wind downstream of the termination shock is 
not well known: the post-shock wind in the head of the neb¬ 
ula, more specifically in the nose region between the termi¬ 
nation shock and the contact discontinuity, is re-accelerated 
to supersonic speeds, while the wi nd crossing the backwar d 
spherical shock is subsonic (see e.g. iBucciantini et al.ll2C)C)^ . 
The average state of the wind is thus not well defined. This 
result seems to be at odds with mass l oading models de- 
veloped in a pure 1-D geometry, (see e.g. ISzego et al.ll200d . 
§2, for a review). More specifically, the solution in the 1-D 
model generally predicts that mass loading in a supersonic 
flow leads to the deceleration and heating of the flow, while 
mass loading in a subsonic flow causes acceleration and cool¬ 
ing. The reason why this behaviour is not observed in the 
present model is due to the fact that the flow in a quasi 1-D 
model can expand in the radial direction while the pressure 
remains constant, and therefore no acceleration in the flow 
direction is possible. On the other hand, if the radial expan¬ 
sion were limited for some reason {e.g. due to hoop stresses 
caused by to a toroidal magnetic field), then an accelera¬ 
tion of the flow would be possible. This specific point will 
be developed in more detail in a subsequent paper. 




Figure 5. The same as Fig. 0 but for scenario B. From top 
to bottom: solutions are calculated for the initial Mach number 
equal to 3, 1 and 0.3, respectively, while the initial wind velocity 
is identical for all panels, he., ui/Vq = 1000. 


Table 2. Summary of the parameter values used for scenarios A 
and B. 


mSM 

A^ion 

Tism 

Vo 

ui Ml 


[cm~®] 


[K] 

[kms~^] 


[ergs-l] 

0.1 

0.9 

lO'^ 

300 

c 1 

1034 
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4 POSSIBLE IMPLICATIONS OF MASS 

LOADING FOR THE BOW SHOCK 

One of the limitations of the present, quasi 1-D model is that 
the bow shock which separates the unshocked ISM from the 
shocked ISM is neglected. It is therefore useful to compare 
the results of the present model with the expected structure 
of the bow shock when mass loading is absent. In Fig. © 
the expansion proHles for scenarios A and B presented in the 
previous section are compared with the profile of an ideal 
bow shock. The bow shock p roHle is calcu lated using the 
thin-shock approximation fsee lWilkinlll996l . Eq.(9)) close to 
the head of the nebula. When the distance from the head 
becomes large, the thin-shock approximation is no longer 
valid, and we therefore replace this solution with the Mach 
cone, where the inclination angle, 6, between the bow shock 
and the pulsar velocity is such that sin0 = I/Mns- For the 
bow shock profile the benchmark values summarised in Ta¬ 
ble [5] are again used. For the mass loaded wind proHles three 
different cases are plotted: scenario A with X/do = 1, and 
scenario B with X/do = 5 and X/do = 30. Fig. (O shows that 
for all mass loaded cases the pulsar wind in the tail of the 
nebula expands beyond the position of the unperturbed bow 
shock, and one would thus expect the geometry of the bow 
shock to be affected. However, while the expanding wind 
profile in scenario A closely follows the proHle of the bow 
shock, the expansion of the wind proHle in scenario B is much 
faster, and is capable of producing the head-shoulder shape 
observed in some Ha bow shock nebulae. It should also be 
noted that the wind profile in scenario A never crosses the 
bow shock when X> do, and one would therefore not expect 
such a configuration to affect the bow shock profile. 

We note that neglecting the primary bow shock {i.e., ne¬ 
glecting the ram pressure of the ISM) leads to an incorrect 
prediction for the wind profile when this profile increases 
beyond the Mach cone. We show in Appendixthat when 
this occurs the formation of a secondary shock in the region 
between the primary bow shock and the contact disconti¬ 
nuity is expected. Such a situation is schematically shown 
in Fig. 0. Behind the secondary shock the pressure will in¬ 
crease, resulting in a bending of the contact discontinuity 
towards the axes of the wind. Consequently the rapid ex¬ 
pansion predicted in scenario B will probably be attenuated 
by the presence of the secondary shock. 

An interesting consequence related to the presence of 
the secondary shock is that the increase in temperature 
will lead to an enhancement of Balmer emission just behind 
this shock. However, the enhancement of the emission will 
strongly depend on the inclination angle of the secondary 
shock: if the inclination angle is close to the Mach cone, the 
temperature and Balmer emission will increase only slightly, 
while for a larger inclination angle the temperature will in¬ 
crease significantly, resulting in a strong enhancement of the 
Balmer emission. 

Although the quasi 1-D model has limitations, it is nev¬ 
ertheless interesting to compare the expansion speed of the 
pulsar wind with the sound speed of the external ISM. This 
expansion speed along the radial coordinate r, as seen by an 
observer co-moving with the external medium, is 

, , dr Vo dA Vo dAdu 

Urix) = — = -;-— = -;-;-— , (40) 

dt du dx 

where we have used A = yrr^ and t = x/Vo- The last equality 


can be used to express Ur in a compact form as a function 
of u only. Using du/dx from Eg. (1331) and calculating dA/du 
from Eq. dSD, leads to 


U 7 . — Vio 


Ai 

47rA^ 


Fin), 


where the dimensionless function F{u) reads 




A' ui(u-Vo)^ 


Ai u^{ui — Vo) 
The value of F at a; = 0 is 




Ml 


Fin,) = ( 1 + 


U\ 


U\ 


(41) 


. (42) 


(43) 


for both scenarios A and B. In the limit Mi <§; 1 and ui 3> 
Vo one has F{u{) = 1, and consequently the value of Ur at 
the origin is 


Url — 




(44) 


By contrast, for very large distances one has lim 2 :_ioo F{u) = 
0. The behavior of the radial expansion speed can be seen in 
Figs.dH) and 0 where the normalised speed Ur{x)/uri = 
F{u) is shown for all plots. With a simple study of the 
function F{u) it is easy to demonstrate that Ur is a mono- 
tonically decreasing function of x for scenario A, while for 
scenario B the function F{u) has a peak where u equals 


_ 5VbMi(2Mi-I-a) 

peak FVoa + ui (2ui + a) ’ 


(45) 


where a is given by Eg. (13811 . In the limit wi ';$> Vq one has 
Mpeak —>■ 5Vb, which corresponds to a value of the expansion 
speed that is equal to 


Ur,peak — Ur(Upeak) — Url A 


16 


Vo 25^5 




1 



3/2 


(46) 

In addition, always in the limit ui ^ Vq, the peak of the 
expansion speed is located at 


^peak 


A 


2 +M2 (7,-1) 


- + 21n 
2 


Ul (2 -)- Ml (7g — 1)) 


8 U 0 


(47) 


Using parameter values that are typical for a bow shock 
PWN leads to Ur,peak ~ (10 — 100)Url and Xpeak ~ (2 — 6 )A. 
Comparison of Fig. © and Fig. ([6ll shows that for A < lOOdo 
the expansion speed in scenario B is larger than the sound 
speed in the ISM when the wind crosses the unperturbed 
bow shock profile. In other words, when the pulsar wind 
expands beyond the unperturbed bow shock, the expansion 
is supersonic and one expects a strong modification of the 
bow shock. By contrast, for A > lOOdo the expansion of the 
wind in the ISM is subsonic, and one therefore expects a less 
pronounced deformation of the bow shock profile. 

From Fig. © it should also be noted that for A < do, the 
velocity Ur can be larger than the sound speed in the wind 
as well as the wind velocity u. When the former condition 
is realised the stationary approach is no longer valid, while 
the second condition leads to a break down of the quasi 1-D 
approximation. As a result our model can no longer be used. 

Based on the above arguments one may speculate that 
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Figure 6. Comparison between three mass loaded wind profiles, 
calculated using the non-relativistic model (solid lines for sce¬ 
nario B and dashed line for scenario A), and the profile of the 
unperturbed bow shock, c alculated usin g the thin shock approx¬ 
imation (dot-dashed line) <IWilkinlll99fih . The results for scenario 
B are shown for two different values of the mass loading distance, 
A = 5do (thin lines) and A = 30do (thick lines), while scenario A 
is plotted only for A = do- All cases have an initial Mach number 
Ml = 1. The shaded region, with radius r = do, shows the mass 
loading region for scenario A. 

when A decreases to a value close to do, the non-stationarity 
of the problem could give rise to a periodic structure of ex¬ 
panding bubbles, similar to those observed in the Guitar 
Nebula. However, addressing this scenario using only ana¬ 
lytical models is a complicated matter, requiring the use of 
full 2-D numerical simulations. 

From the present investigation it is difficult to predict 
which of the scenarios, A or B, is the more realistic one. 
Based on the 1-D approach, one can state that when the 
expansion occurs inside the bow shock, the behaviour is most 
likely well described by scenario A. In this scenario neutrals 
repeatedly undergo charge exchange in the shocked ISM, 
acquiring the same bulk speed, while also flowing parallel to 
the contact discontinuity between the shocked ISM and the 
relativistic wind. On the other hand, when the wind is close 
to the position of the unperturbed bow shock (or crosses 
it) the effective distance between the new bow shock and 
the contact discontinuity could be small enough to allow a 
relevant fraction of neutral particles to penetrate into the 
wind. If this is indeed the case, the expansion speed should 
increase notably, reaching values close to the value predicted 
by scenario B. It should also be noted that far from the 
head of the nebula the neutral fraction of the ISM should 
be larger as the ionization due to the UV radiation emitted 
by the nebula is less effective, and that this should produce 
a faster transition to scenario B. 


5 HYDRODYNAMICAL MODEL FOR 
RELATIVISTIC WIND 

In this section we present the solution of a mass loaded wind 
accounting for the fact that the pulsar wind is relativistic, 
hence we will use a relativistic expression for the enthalpy. 
This is the only difference with respect to the non-relativistic 
treatment presented in §(3] We will show that the solution is 



Figure 7. Sketch of the generation of a secondary shock inside 
the first bow shock due to the mass-loading induced expansion of 
the wind. The shaded region behind the secondary shock shows 
where one would expect an enhancement of the Ha emission. 

very similar to the non-relativistic one, apart from a different 
value of the transition length scale A. 

The relativistic equations for the conservation of par¬ 
ticle number, energy and momentum, written in the rest 
frame of the neutron star and for a 1-D system, are 


dx [up^euA] = hA ', 

(48) 

dx [w7™uA] = qc^oA' , 

(49) 

A] - 1 - c?AdxP = qc?^oA'Vo ■ 

(50) 


Here ne,p is the numerical density of electrons (protons), w 
is the total wind enthalpy and 7 ^, and 70 are the Lorentz 
factors of the wind and the neutrals, respectively. The quan¬ 
tities q and h are given by Eqs. dlTj and (I28II . respectively. 
We notice that in a relativistic treatment it is more conve¬ 
nient to start with the conservation of the particle number, 
Ea. (l48ll . rather than with the conservation of mass. For the 
reader’s convenience, the derivation of Eqs. (gHJ-dSoli is re¬ 
ported in Appendix O 

As was the case for the non-relativistic model, we again 
neglect the free wind region and the termination shock. We 
consider only the pulsar wind after the termination shock 
where it becomes marginally relativistic with a bulk Lorentz 
factor 7 ™ « 1. Furthermore, the neutrals are non-relativistic, 
hence 70 ~ 1 , while the steady-state assumption requires 


the pressure to be constant everywhere, i.e., P = Pq and 
dxP ~ 0. Using these assumptions, it is thus possible to 
simplify the system (Iig1l-(I5U1| as follows: 

doo [peuA] = nm^A' « q{m^/mp)A', (51) 

9a: [ppuA] = firripA' Ri qA' , (52) 

dx [umA] = qp A' , (53) 

dx [uiu^A] = qp A'Vo ■ (54) 


Note that in Eas. (l51ll and (I52II we neglect the contribution 
of the electrons to mass loading, and that the approximation 
q = h(me + mp) ~ hmp is used. In order to close the system 
(ETJ-dSlll, an expression for the enthalpy is required. It is 
generally believed that the pulsar wind predominantly con¬ 
sists of electron-position pairs with highly relativistic tem¬ 
peratures, and we assume that when mass loading occurs 
there is no energy transfer between electron and protons. In 
other words, electrons and protons do not reach a thermal 
equilibrium, but evolve independently with different temper¬ 
atures. This implies that the electron gas is always highly 


© 0000 RAS, MNRAS 000, 000-000 






12 G. Morlino, M. Lyutikov and M. J. Vorster 


relativistic, hence the rest mass contribution to enthalpy 
can be neglected, i.e., We = f-e + Pe = 4Pe = 4Po- On the 
other hand, protons are non-relativistic, hence their thermal 
energy is always negligible with respect to their rest mass 
energy, and their enthalpy is Wp = Pp(? ■ The total enthalpy 
of the wind is thus 

W = Wp We = PpC + 4Po • (55) 

Using these simplifications allows one to obtain an analyti¬ 
cal solution for the wind dynamics by following a procedure 
similar to the one outlined for the non-relativistic case. Com¬ 
bining Ea. (l53ll and (1541) furnishes a first constant of motion: 

uAw{u — Vb) = uA{ppC? + 4Po)(m — Vb) = const. (56) 

Similarly, combining Eg. (1521) and Eg. (1531) furnishes a second 
constant of motion: 

uA(w — PpC^) = 4uAPo = const. (57) 


Evaluating the constants at the initial position x = xi, 
where u = ui and A = Ai, Eg. (1571) gives the solution for 
the cross section ^4 as a function of the velocity: 


A{u) = uiAiju , 


(58) 


while dividing Eg. (1561) by Eg. (1571) gives the solution for the 
proton density as a function of u: 


Pp{u) 


4Po Ui — u 
(? u — Vo ' 


(59) 


The electron density can be obtained in a similar way using 
Eo. (|51|) rather than Eo. (l52l) . giving the following result 


Pe(u) = pe,0 + 


me 4Po ui — u 
mp u — Vo ’ 


(60) 


where pe,o is the initial electron density. One noteworthy 
point is that the solutions (1551) . (15^ and dnni) do not depend 
on the specific assumption regarding the injection cross sec¬ 
tion A' (this result is also identical to the non-relativistic 
case). The solution for the velocity u{x) can be obtained by 
deriving Eq. (I5il) by parts, and using Eq. (1551) . This leads to 
the expression 


du _ 2 M — Vb 

dx A wu 


(61) 


As was done for the non-relativistic case, we again dis¬ 
tinguish between scenario A {A' = Ai) and scenario B 
{A' = A). For scenario B the integration of Eg. (1611) leads 
to the following implicit solution: 


X 

Arel 


u 

u —Vo 


m 

ui — Vo 


4-In 


m — Vo 
M — Vb 


(62) 


where we have introduced the relativistic length scale 

4Po(ui —Vb) 4Po Ml , . 

Arel = -5- — - p -r- • (63) 

qc^ Pnc^ UphCTphC 

The second equality results from assuming Vb ^ ui ~ c and 
using the photo-ionization rate riphitphC in the calculation of 
the mass loading rate, i.e., q = mpUNriph^phC. Comparison 
of Eg. (1631) with the definition of the mass loading length 
scale for the non-relativistic case, Eq. (O. shows that they 
are identical, with the exception that the quantity 4Po/c^ 
replaces the electron mass density. In the steady-state model 
the internal pressure of the wind is equal to the external 
pressure, hence 4Po corresponds to the specific enthalpy of 


the electron-positron wind plasma. Using Eq. Hi) from 1)2.31 
allows one to find a numerical estimate for AreP 

Arel « 1.2 • 10“ T4^ ^ n^;_ 4 .E(r)"^Cm, (64) 

where Po = niswiKsT and T = 10'*r4 K have been used. 
Choosing realistic values for the parameters associated with 
bow shock nebulae emitting Ha shows that Arei can be as 
large as 10^® cm, but we stress that this value can vary by 
orders of magnitude, essentially due to the fact that the val¬ 
ues of the neutral density inside the wind and the luminosity 
of the PWN tail are difficult to estimate. Nevertheless, for 
any scenario Arel represents the lower limit for the typical 
expansion length scale. For example, in scenario A a larger 
expansion length scale is predicted. 

Substituting A' = A\ in Eg. (1611) leads to the following 
equation 


du _ qc^ (m — Vb)^ 
dx 4Po (ui — Vo)ui ’ 


which, once integrated by parts, gives the solution 

X Ul Ul 

Arel u — Vo Ml — Vb 

In this case the typical expansion length scale is 


^ Ml ^ 4 Po(mi - Vb)Mi 

Arel,A — -r r Arel — orr 

Vo qc^Vo 


(65) 


( 66 ) 


(67) 


As Mi/Vb can be as large as 10®, the mass loading length 
scale in scenario A can reach 10*^® cm. In a realistic case one 
expects the transition to occur between Arei and Arei.A- 

Eg. (1641) predicts that a density of neutrals as small as 
10“'* cm“® is sufficient to have Arei comparable to the size 
of the nebula. Thus, for a relativistic pair-plasma wind, a 
relatively small amount of neutrals can strongly affect the 
tail flow. 

Fig-® shows the solutions for u and A corresponding 
to scenarios A and B, where the parameter values given in 
Table[2]have been used. Comparison of Fig.® with Figs.® 
and ® shows that the relativistic and non-relativistic solu¬ 
tions are very similar. 

As was done for the non-relativistic case, the profile of 
a relativistic wind undergoing mass loading is compared to 
the typical profile of the unperturbed bow shock. Fig. ® 
is analogous to Fig. ®, but the wind profile is now calcu¬ 
lated using the relativistic expression for scenarios A and 
B obtained in Eg. (1581) . One can see that the relativistic and 
non-relativistic model give very similar profiles, with the im¬ 
portant difference that A is replaced by Arei- Additionally, a 
similar radial expansion is predicted for the corresponding 
scenarios of the non-relativistic and relativistic cases. Start¬ 
ing from Eg. (1401) . using du/dx from Eg. (1651) . and calculating 
dA/dx from Eg. (1581) . leads to the following expression for the 
expansion speed: 


Ur — Vb ^ 


/ 

M^^^ (m - Vb)^ 

\J 47rA2^i 

4(5/2 



( 68 ) 


where the upper and lower signs refer to scenarios A and 
B, respectively. Fig. ® also shows the normalised velocity 
Mr/Mri, where 


Mt-i = Vb 



do 

2Ai.el 


4b, 


(69) 
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Figure 8. Structure of a pulsar wind when mass loading occurs, 
assuming that the wind consists of a hot relativistic electron- 
positron pair plasma. Top and bottom panels report results for 
scenario A and B, respectively, using the same values summarised 
in Table [21 The various lines represent the wind speed divided by 
Vo (solid line), wind cross section (dot-dashed line) and expansion 
speed (dotted line) both normalised to their initial values. The 
sound speed of the ISM, Cso, also normalised to Vq? is shown for 
comparison (thin dot-dashed line). 



Figure 9. Same as Fig.©, but for a relativistic model. Notice 
that the dot-dotted line is the same as in Fig. © and represents 
the bow shock profile in the thin shock approximation for non- 
relativistic wind. 


is the expansion speed at a: = 0, with this speed being 
the same for both scenarios A and B. Similar to the non- 
relativistic case, Ur(x) is also a monotonically decreasing 
function of x for scenario A, while for scenario B it has a 
peak at M = Mpeak = 5Vb, corresponding to 


ttr.peak — ‘^r (’^peak) — Url i 


16 


Vo 25V5 ' 


(70) 


The peak is located at the position Xpeak, which, in the limit 
ui Vb, is equal to 


^peak -^rel 




(71) 


It is interesting to note that this result is identical to the 
non-relativistic case for the limit Mi —^ 0. In other words, 
all discussions presented for the non-relativistic solutions at 
the end of (©also apply to the relativistic solutions. In par¬ 
ticular, when Arei ^ do the expansion velocity can become 
larger than the wind velocity, thereby causing the quasi 1-D 
approximation to break down. 


6 A VISUAL FIT OF BOW SHOCK NEBULA 

PSR J0T42-2822 

In this section we provide an example of a comparison be¬ 
tween the predicted profile of the wind due to mass loading 
and the observed shape of an Ha bow shock. Among the 
pulsars showing an anomalous Ha bow shock, we chose PSR 
J0742-2822 for two reasons: not only has several quantities 
been measured with sufficient accuracy, but the pulsar is 
believed to move almost perpendicular to the line of sight, 
a fact which, to some extent, simplifies the comparison be¬ 
tween the model prediction and observation. 

We use the pa ra meter values reported by 
iBrownsberger fc Romanll ll20l4l . The pulsar luminos¬ 
ity is Cm = 1.9 • lO^^ergs”^, while the measured proper 
motion is 29.0 mas yr~^ . The most recent distance estimate 
gives 2.1 ± 0.5 kpc (|ja.nssen fc Stappersl l2006li. Using 
a value of d = 2 kpc, Brownsberger fc Romanil 


estimated nisM = 0.28 cm“'^ for the total ISM density and 
V± = 275kms“^ for the projected component of the pulsar 
velocity. We assume that the velocity of the pulsar is almost 
perpendicular to the line of sight, hence Vo ~ V±. 

Using the above values, the stand-off distance is esti¬ 
mated to be do = 3.8 ■ 10^® cm. Note that this value is 
compatible with the measured distance between the pulsar 
and the apex, measured as 1.4", corresponding to a phys¬ 
ical distance of 4.2 • 10'^®d/(2kpc) cm. Using this value of 
do, the profile of an ideal bow shock (ie., no mass load¬ 
ing) is calculated according to the procedure outlined at the 
beginning of © The resulting solution is shown in Fig. [10] 
(dot-dashed line), while the small dot indicates the position 
of the pulsar. The unperturbed bow shock profile fits the 
head of the nebula reasonably well, but fails to reproduce 
the “fan” structure emerging behind the bow shock head. 

For the mass loaded wind, profiles are calculated using 
the relativistic model developed in © with Fig. [10] showing 
the results for both scenarios A (dashed line) and B (solid 
line). In order to ensure that the wind profile in scenario 
B crosses the bow shock at the precise location where the 
“fan” structure emerges, the profiles are calculated using 
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the value Arei = 2do- Note that the wind profiles start at the 
position dback = 4do behind the pulsar, which corresponds 
to the location of the bac kward termination shock according 
to numerical simulations dBucciantini et al]l2005ll . As previ¬ 
ously discussed in ^and ^ scenario A cannot account for 
a rapid expanding structure emerging from the unperturbed 
bow shock as the resulting profile is very smooth, while sce¬ 
nario B does predict such a structure. On the other hand, 
from Fig.[Tn]it can be seen that the expansion in scenario B 
is larger than the one detected in Ha. 

Our aim is not to exactly reproduce the shape of the 
bow shock as our model only predicts the shape of the rel¬ 
ativistic tail wind. The more important point that we want 
to emphasise is that the location where the “fan” structure 
emerges is compatible with our predictions. In particular, 
one can ask whether the ratio Arei/do = 2 is compatible 
with observations. This can be checked by estimating the 
neutral density inside the wind. _ 

Analysing the Ha flux. lBrownsberger fc Romanil il2014l l 
found that the neutral fraction of the ISM is comparable to 
~ 1. However, as the uncertainty is quite large, we use a 
fiducial value of 0.5. It was shown in that the colli- 
sional ionization in the shocked ISM is not effective when 
Vb < 300kms“^, therefore one need only take into account 
the effect of photo-ionization. As was shown in >12.31 the UV 
flux from the nebula can be estimated from the X-ray lu¬ 
minosity: the measured upper li mit for the X-ray flux is 
Fx < 2 • erg cm“^ s“^ llAbdo et al.l 1201 j f. which 

translates into an upper limit for the X-ray luminosity of 
Cx < 9.6 • 10®'’(d/2kpc)^ ergs“^. For the non-thermal spec¬ 
trum we assume a power law with an index F = 2. Further¬ 
more, using Eq. m allows one estimate the neutral density 
at the location of the backward termination shock, dback- 
The value of this density averaged over the cross section of 
the wind is found to be ~ 0.13 nisM = 0.065 cm“®. Finally, 
using this value in Eg. (1 641) and dividing the result by Eq.([T]), 
one finds that Arei/do > 0-03R/{l). A value of Arei/do = 2 
thus implies R/{1) < 60, which is fully compatible with the 
wind geometry. 


7 MASS LOADING IN A MAGNETISED WIND 

In this paper we have neglected the role of the magnetic 
field inside the wind of the nebula. Nevertheless, using the 
conservation of magnetic flux, we now show that general 
conclusions can be drawn from the knowledge of the wind 
dynamics only. We consider two different configurations: a 
purely poloidal and a purely toroidal magnetic field, with 
both cases assuming that the initial magnetic field is dy¬ 
namically unimportant. In the poloidal case the flux conser¬ 
vation is applied through the cross section A, which leads 
to BxA = const, where is the average poloidal magnetic 
field strength. Using the solution for the cross section from 
Eg. (1581) . one finds that Bx oc u. As the fluid speed decreases, 
the poloidal magnetic field strength will also decreases, and 
one would not expect the wind dynamics to be affected. 

For the toroidal configuration the opposite situation oc¬ 
curs. The conservation of magnetic flux is applied through 
a poloidal surface (i.e., a surface parallel to the x direc¬ 
tion) which implies B^uy/A = const, where B^ is the av¬ 
erage poloidal magnetic field inside the wind. Again using 



Figure 10. Comparison of the observed Ha bow shock produced 
by the PSR J0742-2822 l lBrownsberger fc Romanill2014h with the 
wind profile obtained using the relativistic calculations for sce¬ 
nario A (dotted line) and scenario B (solid line). The dot-dashed 
line shows the predicted position of the unperturbed bow shock. 
The small dot in the head of the nebula indicates the position of 
the pulsar. The profile of the wind starts at the location of the 
backward termination shock, which is Ado behind the pulsar. 

Eg. (1581) . one Hnds B^ oc Consequently the strength 

of the toroidal component increases along x. As the toroidal 
magnetic held exerts a hoop stress on the wind oc B^, the 
magnetic pressure can easily become comparable to the in¬ 
ternal pressure, thereby reducing the expansion of the wind. 

These conclusions have implications for the non-thermal 
emission. The synchrotron emissivity is proportional to 
jsyn oc neB^, where rie is the density of relativistic elec¬ 
trons and B is the average held. The density of relativistic 
electrons is constant along x, and as a result the expansion 
induced by mass loading will enhance jsyn if the magnetic 
held is mainly toroidal. Conversely, for a poloidal conhgura- 
tion at the location of expansion the synchrotron emission 
should decrease. 

Remarkably, this prediction is compatible with the ob¬ 
servations of at least one object. The bow shock nebula pow¬ 
ered by the pulsar J1509-5850 has been detected both in X- 
rays an d radio, and shows a n anti-correlation in th ese two 
bands liHui fc Becked l2007l : iKargaltsev et al.l l2008lj : while 
the head of the nebula is mainly bright in X-rays, the bulk 
of the radio emission is observed from the tail of the wind at 
a distance roughly 4’ far away from the pulsar, at a location 
where the X-ray emission becomes negligible. Moreover, ra¬ 
dio polarimetry measurement s show that the magnetic field 
in the tail is mainly toroidal (iNg et al.ll201Cll l. 

In our model these observations can be explained by 
assuming that the enhancement of the radio emission oc¬ 
curs where the tail expands due to mass loading, while the 
decrease in X-ray emission should be the consequence of ra¬ 
diative cooling, which is efficient for the X-ray emitting elec¬ 
trons but is negligible for the radio emitting electrons. This 
picture is conhrmed by the fact that J1509-5850 has also 
been observed in Ha, showing an anomalous bow shock ex¬ 
pansion right at the location where the enhancement of the 
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radio emission begins. A further intriguing detail is that the 
comparison between radio and X-ray data suggests a signifi¬ 
cant deceleration of the flow when moving downstr eam, with 
the s peed decreasing by 1-2 orders of magnitude dNg et al.l 
l2010lf . From our model it follows that this decreases can be 
the result of mass loading. However, a more quantitative 
model is needed to better compare this picture with obser¬ 
vations. 

Noticeably, an opposite situation has been observed 
for the Mouse nebula (J1747-2958), where radio polarime- 
try shows a poloidal magnetic field structure along the tail 
dGaensler et al.ll2004h . In this case X-ray and radio emis¬ 
sion are both peaked in the head of the nebula and decrease 
smoothly along the tail. The effect of mass loading for the 
Mouse nebula is probably negligible, in fact the Ha emission 
is not detected, and this is probably due to its exceptional 
luminosity able to photo-ionise the ISM around it. The X- 
ray luminosity measured b y Chandra is Cx = 5 • 10®'* erg 
s“^ for a distance of 5 kpc (iGaensler et al.|[2004ll . while the 
pulsar spin down luminosity is Cw ~ 2.5 x 10^^ e r g s~^ , 
giving an X-ray efficiency = 0.02. IGaensler et al.l (l2004l l 
also measured the photon index as 1.8 < Fx < 2.5, and pro¬ 
vide an estimate of the ISM density of ~ 0.3 cm“^. Using 
these values we can estimate the upper limit for the X-ray 
efficiency from Eq. (HI to be 1.2-10 ■*. Because r]x is much 
larger than this upper limit, we infer that the vast majority 
of H atoms have been ionized before to reach the bow-shock. 


8 DISCUSSION AND CONCLUSIONS 

In this work we have shown that the structure of a bow shock 
nebula produced by a neutron star propagating through the 
ISM can be significantly affected by the presence of neutral 
hydrogen in the ISM. In many cases a non-negligible fraction 
of neutral atoms can penetrate into the relativistic wind 
where they will be ionized by UV photons emitted by the 
nebula. Once Ionized, the new protons and electrons interact 
with the wind, leading to a net mass loading of the wind in 
the tail region of the nebula. More specifically, this mass 
loading is important in the shocked part of the tail wind, 
while it is most likely negligible in the region of the free, 
unshocked wind. 

To investigate the effect of mass loading, we have de¬ 
veloped a steady-state hydrodynamic model using a quasi 
1-D approximation, with the study focusing on two specific 
scenarios. In the first a wind consisting of a non-relativistic 
plasma with an adiabatic index 7 g = 5/3 (see 331 was in¬ 
vestigated, while the second situation investigated a wind 
that consists of a hot relativistic electron-positron plasma, 
as is the case for a pulsar wind (see 33- Remarkably, both 
situations show the same qualitative behavior: the loaded 
mass decelerates the wind and produces a transverse expan¬ 
sion of the tail. The typical expansion factor of the tail cross 
section is Aoo/Ai = ui/Vb, where ui ~ c is the initial ve¬ 
locity of the pulsar wind after the termination shock, and 
Vo is the speed of the neutron star measured in the rest 
frame of the ISM, which is of the order of few hundreds 
kms“^. The model therefore predicts Aoo/Ai « 1000 (for a 
non-relativistic wind there is a correction due to the initial 
Mach number of the wind - compare Ea. (l39ll with Eg. (1581) 1. 
The main difference between the two situations is the dis¬ 


tance where this expansion occurs, being determined by the 
mass-loading length scale, Aml. We further showed that this 
latter length scale is determined by the distance where the 
enthalpy of the loaded mass is comparable to the enthalpy 
of the wind (compare Eq. (fTfl) vs. Eq. (l63tl. 

Considering a relativistic wind consisting of an electron- 
positron plasma, and using parameters values observed for 
pulsar bow shock nebulae, we showed that the mass loading 
effect could be responsible for the head-shoulder shape ob¬ 
served in many bow shock nebulae. Unfortunately it is not 
easy to make a quantitative prediction for specific objects as 
the amount of mass loading depends on two important pa¬ 
rameters, specifically the density of neutrals inside the wind 
and the density of UV photons, both of which are difficult 
to estimate. 

We showed that a relatively small density of neutrals 
inside the wind (as small as 10“‘^cm“®) is sufficient to af¬ 
fect the wind, decelerating the tail flow and producing a 
fast expansion in the transverse direction. For comparison 
we remember that the typical number density of neutral 
hydrogen in the warm interstellar medium (where Ha bow 
shock nebulae are though to propagate) is ~ 0.05 — 0.5 cm“® 
djean et al.ll2009ll . In order for the mass loading effect to be 
negligible, either the ISM should be completely ionized, or 
the photo-ionization due to the nebula should be so effective 
that all neutrals are ionized before they reach the termina¬ 
tion shock in the tail. Neither of these two possibilities can 
be applied to Ha bow shock nebulae as the presence of Ha 
lines is an indication that the photo-ionization is incapable 
of fully ionizing the ISM. If this were the case, it would not 
be possible to observe the Ha emission in the first place. 

We also made a visual comparison of the wind profile 
predicted from our model with the Ha bow shock observed 
around PSR J0742-2822. We showed that a density of neu¬ 
trals in the wind comparable to ~ 0.06 cm“® is required to 
explain the presence and the location of the “fan” structure 
observed behind the head of the nebula. Such a density is 
compatible with the estimated neutral density of the local 
ISM, taking into account the photo-ionization resulting from 
the UV radiation emitted by the nebula, which, in turn, was 
estimated from the upper limit of the X-ray flux. 

Finally, even though the magnetic field is neglected in 
the present model, we discussed the qualitative behaviour of 
a magnetised wind flow that is being loaded with mass. The 
poloidal component of the field decreases with distances, 
thereby becoming dynamically unimportant. On the other 
hand, the toroidal component, if present, will be amplified 
due to the compression of the wind flow. Such an ampli¬ 
fication will have two important consequences: the first is 
to limit the transverse expansion of the wind due to the 
increase of the magnetic hoop stresses; the second is to pro¬ 
duce an enhancement of the synchrotron emission in the 
location where the expansion occurs. Remarkably, this pre¬ 
diction is confirmed by the observation of the bow shock 
nebula associated with the pulsar J1509-5850, where an en¬ 
hancement of the radio emission is observed far from the 
head of the nebula, and at the same location where the Ha 
emission shows an anomalous expansion of the bow-shock. 

The quasi 1-D approximation presented in this work has 
the advantage of having a clear and simple analytical solu¬ 
tion, thereby making it very usable. However, it also has 
several limitations. Firstly, we neglected the presence of any 
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internal structure, as well as the bow shock. Moreover, the 
quasi 1-D approximation is based on the assumption that 
the transverse expansion velocity is much smaller than the 
wind velocity, a condition that can be violated when Aml 
is of the order of, or smaller than the typical stand-off dis¬ 
tance of the nebula. In this case the quasi 1-D approximation 
breaks down. Additionally, the transverse expansion speed 
can be faster than the sound speed of the wind. When this 
occurs the steady-state condition is violated and the solution 
presented is no longer applicable. It is therefore necessary to 
implement a time-dependent solution. This last situation is 
especially interesting as it points toward to possibility of 
having a tail wind with periodic expanding bubbles simi¬ 
lar to what is observed in the Guitar nebula. Hopefully all 
this limitations can be overcome using time-dependent, 2- 
D simulations, which will form the topic of research of a 
forthcoming paper. 


^(a;), relating the new variable ^ to the physical distance x. 
This can be easily done using Eq. (IMl) along with the equa¬ 
tion for the evolution of the neutral density, which, for the 
scenario A, is: 

dUN . _ , . 

Vo-^ =-n =-njvnphCTphC, (A3) 

leading to a simple exponential solution 

riNix) = , (A4) 

where Aph is defined in Eq. (II6J. Finally, substituting 
Ea. (IA4ll into Ea. IAlll . and integrating one obtains 

C(») = Aph . (A5) 

It is easy to check that when x -C Aph one has ^ = x, and 
that the solutions presented in ^and !j5]are recovered. 
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APPENDIX A: DEPLETION OF NEUTRALS 
INSIDE THE WIND 

For the solutions presented in the paper we have assumed 
that the density of neutrals inside the wind remains con¬ 
stant, i.e., that the number of ionized atoms is negligible 
with respect to the total number of neutrals. This is true 
only on a length scale much smaller than the ionization 
length scale, Aph. Here we show how to modify the model, 
should this assumption no longer be valid, by taking into 
account the fact that the density of neutrals is no longer 
constant, but decreases as a function of distance. 

When the depletion of neutrals becomes important, sce¬ 
nario B (where neutrals can also penetrate into the wind 
region through the side of the contact discontinuity) can no 
longer be described using a 1-D model as the neutral density 
inside the wind will depend on both x and r. The discussion 
will therefore be limited to scenario A, where the neutrals 
penetrate only through the head of the nebula. 

It is useful to introduce a new spatial variable, defined 
as 

= nN{x)dx, (Al) 

where un = nN{x)/nNi is the numerical density of neu¬ 
trals normalised to its initial value. Substituting ^ for nN{x), 
Egs. (1241) - (1261) (or Eas. (l48l) - (I50I) for the relativistic case) can 
be rewritten by replacing ^ d(^, and using 

hi = nivinp/io-phc (A2) 

rather than h on the right-hand side. This new set of equa¬ 
tions is formally identical to Eqs. where a constant 

density is assumed, with the only exception that x is re¬ 
placed by The solutions presented in !j3]and ^therefore 
remain the same, with the exception that the dependence on 
X is substituted by the dependence on To find the full new 
solutions therefore only requires that one finds the function 


APPENDIX B: FORMATION OF SECONDARY 
SHOCKS 

One of the main results of this study is that the mass loading 
of a pulsar wind leads to a sideways expansion of the tail. 
As the contact discontinuity between the tail flow of a PWN 
and the ISM acts as impenetrable wall for the ISM flow, a 
secondary shock can be created. As we demonstrate below, 
this secondary shock will appear when the contact disconti¬ 
nuity makes an angle with the flow that is larger than the 
Mach cone. 

As a model problem, consider a flow with velocity, V, 
and internal sound speed, Cs, bounded by a wall with the 
parabolic profile, rg = XofL. A shock will appear whe n the 
characteristics first intersect llLandau fc Lifshitall959| ). The 
characteristics originating at point {a;o,ro} are 

r{t) = ro -f Cst, 
x{t) = xo + (cs -I- V)t. 

Eliminating t leads to 

r = xl/L-\ -2^(a;-a:o). (B2) 

Cs + V 

The first intersection occurs where dr/dxg = 0: 

^ 1 

2(1-bM) ’ (B3) 

M = V/cs. 

At this point the angle that the wall makes with the flow is 

tan6» = 1/(1-b M). (B4) 

Thus, for a highly supersonic flow a shock forms when the 
angle of attack becomes larger than the Mach angle. 


APPENDIX C: HYDRODYNAMICAL MODEL 
FOR RELATIVISTIC WIND 

To derive Eqs. (|48l)-(|50l) we start by writing a relativistic 
formulation of the mass-loading problem, and then we spe¬ 
cialise our equation s to describe a typ ical pulsar wind. It 
has b een shown by iKomissaro 3 (Il994) (see also lLvutiko3 
l2003l) that the covariant relativistic equation for a perfect 
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fluid with the inclusion of mass loading can be written in 
the following form: 


dx'' 


qcT 


(Cl) 


where 


r"" = wu'n" + Pg"" (C2) 

is the energy momentum tensor, w = e + P is the total 
enthalpy, and = 7 u,(l,u/c) is the four-velocity of the 
plasma, with 7 ™ denoting the Lorentz factor of the wind. 
The quantity qcr^ represents the mass loading term, where 
= 70 ( 1 , Vo/c) is the four-velocity of the neutrals moving 
with a Lorentz factor 70 = (1 — Vo/c)~^^^, and 

q = ^ Ui^irni (C3) 


is the mass injected per unit time per unit volume calculated 
in the frame moving with velocity Vo, i.e., the rest frame of 
the neutrals. It is assumed that the incoming neutrals have 
the same temperature as the ISM (« 10^ K), and they can 
thus be considered as being cold with respect to the wind in 
the nebula, i.e, 'yi = 1. As the rate of injected electrons and 
protons is the same, = hp = h, the equations describing 
the evolution of the electron and the proton densities are 
similar, i.e., 

^ [ne,pcu] = h , (C4) 

where n^^p is the numerical density of electrons (protons), 
with an expression for h given by Ea. (l28ll . As was done 
for the non-relativistic case, we assume that both un and 
Uph are constant along x. A method to obtain the solution 
when the depletion of neutrals inside the wind is taken into 
account is described in Appendix m 

The next step is to rewrite the conservation equations 
m and (IC4I) for a 1-D system, integrating ove r a small vol¬ 
ume with a cross section A and a length dx fsee lKomissarovI 
Il994h . Additionally assuming that the system is in a steady- 
state, the conservation equations for the particles’ number, 
energy and momentum in the x direction become: 

dx [up^euA] = hA', (C5) 

dx [w 7 ™mA] = qP^oA ', (C 6 ) 

dx [wu^A] -I- PAdxP = qP^oA'Vo . (C7) 

Note that the non-relativistic equations (l24l)-([26l) can be re¬ 
covered from Eqs. (IC5l)-(IC7ll using the first order approxi¬ 
mation for the Lorentz factors, i.e., 7 ™ Ri 1 -|- tp je? and 
70 ~ 1 + Vf/c\ 
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